####################################
# Causal mechanisms RDROBUST
####################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/05_causal_mech_crime.txt")

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)

################
# Prepare data 
################

# read data
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

# check transfers
tapply(d$total_tranfers_deflactor,d$year,mean)
mean(d$total_tranfers_deflactor, na.rm = T)
round(mean(d$total_tranfers_deflactor, na.rm = T)/mean(d$amount_projects_deflactor, na.rm = T),2)

###################################
# Effect of alignment on tranfers
###################################

# bandwidth
h = rdbwselect(d$total_tranfers_deflactor_per_s,d$margin,cluster = d$cluster)$bws[1]
round(h,3)

# observations 
sum(abs(d$margin)<=h)

# Total transfers 
summary(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,cluster = d$cluster,all=TRUE))

# Plot
cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figure3.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$total_tranfers_deflactor_per_s, x = d$margin,  h= 0.131, nbins = 100, subset = -0.131 <= d$margin & d$margin <= 0.131, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black",  y.lim = c(-0.5, 2), title = "                                       Robust CI: [0.051 , 0.683]",
        y.label = "Discretionary tranfers", x.label = "Margin of victory"))
dev.off()

############################
# Electricity
############################

# consumo electricidad servicios comunitarios
summary(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,cluster = d$cluster, all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figure4.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$cons_elec_serv_com_def_per_s, x = d$margin,  h= 0.088, nbins = 100, subset = -0.088 <= d$margin & d$margin <= 0.088, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-1, 1), title = "                                       Robust CI: [-0.028 , 0.321]",
        y.label = "Municipal electricity spending on services provided to the community", x.label = "Margin of victory"))
dev.off()

###################################
# More mechanisms 
###################################

# Gastos Municipales Área de Gestión: Programas Sociales 
summary(rdrobust(d$social_programs_per_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA4b.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$social_programs_per_s, x = d$margin,  h= 0.199, nbins = 100, subset = -0.199 <= d$margin & d$margin <= 0.199, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-0.5, 2), title = "                                       Robust CI: [-0.454 , 0.193]",
        y.label = "Social programs", x.label = "Margin of victory"))
dev.off()

# Personas Inscritas en la Municipalidad en Busca de un Empleo 
summary(rdrobust(d$employment_per_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA4a.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$employment_per_s, x = d$margin,  h= 0.038, nbins = 100, subset = -0.038 <= d$margin & d$margin <= 0.038, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-0.1, 0.1), title = "                                       Robust CI: [-0.031 , 0.005]",
        y.label = "Employment", x.label = "Margin of victory"))
dev.off()

#####################
# Table 5
#####################

summary(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,cluster = d$cluster,all=TRUE))
pe1 = as.numeric(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv1 = as.numeric(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci1 = as.numeric(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss1 = as.numeric(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess1 = as.numeric(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw1 = as.numeric(rdrobust(d$total_tranfers_deflactor_per_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

pe = c(pe1)
pv = c(pv1)
ci_lower = c(ci1[1])
ci_upper = c(ci1[2])
oss = c(oss1)
ess = c(ess1)
bw = c(bw1)
variable = c("Discretionary funds")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/table5.html")

#####################
# Table 6
#####################

summary(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,cluster = d$cluster, all=TRUE))
pe1 = as.numeric(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv1 = as.numeric(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci1 = as.numeric(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss1 = as.numeric(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess1 = as.numeric(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw1 = as.numeric(rdrobust(d$cons_elec_serv_com_def_per_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

pe = c(pe1)
pv = c(pv1)
ci_lower = c(ci1[1])
ci_upper = c(ci1[2])
oss = c(oss1)
ess = c(ess1)
bw = c(bw1)
variable = c("Municipal electricity spending")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/table6.html")

#############################
# Table A15 (last two rows)
#############################

summary(rdrobust(d$employment_per_s,d$margin,cluster = d$cluster,all=TRUE))
pe1 = as.numeric(rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv1 = as.numeric(rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci1 = as.numeric(rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss1 = as.numeric(rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess1 = as.numeric(rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw1 = as.numeric(rdrobust(d$employment_per_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

summary(rdrobust(d$social_programs_per_s,d$margin,cluster = d$cluster,all=TRUE))
pe2 = as.numeric(rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv2 = as.numeric(rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci2 = as.numeric(rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss2 = as.numeric(rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess2 = as.numeric(rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw2 = as.numeric(rdrobust(d$social_programs_per_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

pe = c(pe1,pe2)
pv = c(pv1,pv2)
ci_lower = c(ci1[1],ci2[1])
ci_upper = c(ci1[2],ci2[2])
oss = c(oss1,oss2)
ess = c(ess1,ess2)
bw = c(bw1,bw2)
variable = c("Employment","Social programs")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/tableA15_part2.html")

sink()
